package scu.maqiang.homogenization;

import scu.maqiang.fes.*;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.*;

public class SOTS2DT361 {
	public SOTS2DT361(FES2T31 fsT3, FES2T61 fsT6) {
		this.fsT3 = fsT3;
		this.fsT6 = fsT6;
	}
	/**
	 * ��������ȴ��������ͷ����ȴ���ϵ����һ�ס����׵������Լ����Ȼ��ȴ���ϵ����
	 * 
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 * @param SOCS
	 */
	public void KxxCellSolution(MatrixFunc kxxFunc, double[][][] param, double[][] FOCS, double[][] D,
			double[][] SOCS) {
		SRMatrix A = new SRMatrix(fsT3.GetNdof());
		fsT3.assembleStiff(kxxFunc, param, BVPType.CONSTITUTE_MATRIX_COMMON, A);
		double[][] RHS = new double[2][fsT3.GetNdof()];
		fsT3.assembleGlobalVector(kxxFunc, param, SOTS2DT361::elementkxxFirstOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
		fsT3.applyBC_MBN(A, Direct.All, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3, 4);
		MVO.fill(FOCS, 0.0);
		DirectSSolver solver = new DirectSSolver(A);
//		NewIterSSolver solver = new NewIterSSolver(A);
		FOCS[0] = solver.CholeskySolve(RHS[0], 1);
		//solver.PLUSolve(RHS[0], FOCS[0], 1);
		FOCS[1] = solver.CholeskySolve(RHS[1], 1);

		computeHomogenizedKxx(kxxFunc, param, FOCS, D);

		RHS = new double[4][fsT3.GetNdof()];
		fsT3.assembleGlobalVector(kxxFunc, param, FOCS, D, SOTS2DT361::elementKxxSecondOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
		MVO.fill(SOCS, 0.0);
		for(int i = 0; i < 4; i++) {
			fsT3.applyBC_MBN(RHS[i], Direct.All, 0.0, 1, 2, 3, 4);
			SOCS[i] = solver.CholeskySolve(RHS[i], 1);
		}
	}

	public double[][] AijFirstOrderCellSolution_DirichletBoundary(MatrixFunc kxxFunc, double[][][] param) {
		double[][] FOCS = new double[2][fsT3.GetNdof()];
		SRMatrix A = new SRMatrix(fsT3.GetNdof());
		fsT3.assembleStiff(kxxFunc, param, BVPType.CONSTITUTE_MATRIX_COMMON, A);
		double[][] RHS = new double[2][fsT3.GetNdof()];
		fsT3.assembleGlobalVector(kxxFunc, param, SOTS2DT361::elementkxxFirstOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
		fsT3.applyBC_MBN(A, Direct.All, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3, 4);
		DirectSSolver solver = new DirectSSolver(A);
		solver.Factorize_Cholesky(1);
		FOCS[0] = solver.Solve_Cholesky(RHS[0], 1);
		FOCS[1] = solver.Solve_Cholesky(RHS[1], 1);
		return FOCS;
	}

	public double[][] AijFirstOrderCellSolution_PeriodicBoundary(MatrixFunc kxxFunc, double[][][] param) {
		double[][] FOCS = new double[2][fsT3.GetNdof()];
		SRMatrix A = new SRMatrix(fsT3.GetNdof());
		fsT3.assembleStiff(kxxFunc, param, BVPType.CONSTITUTE_MATRIX_COMMON, A);
		double[][] RHS = new double[2][fsT3.GetNdof()];
		fsT3.assembleGlobalVector(kxxFunc, param, SOTS2DT361::elementkxxFirstOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
		fsT3.applyBC_Periodic(A, RHS[0], Direct.Y, 1, 3);
		fsT3.applyBC_Periodic(A, RHS[0], Direct.X, 2, 4);
		fsT3.applyBC_Periodic(RHS[1], Direct.Y, 1, 3);
		fsT3.applyBC_Periodic(RHS[1], Direct.X, 2, 4);
		NewIterSSolver solver = new NewIterSSolver(A);
		solver.PCGSSOR(RHS[0], FOCS[0], 1.5, 1);
		solver.PCGSSOR(RHS[1], FOCS[1], 1.5, 1);
		return FOCS;
	}

	public double[][] AijFirstOrderCellSolution_PeriodicDirichletBoundary(MatrixFunc kxxFunc, double[][][] param) {
		double[][] FOCS = new double[30][fsT3.GetNdof()];
		SRMatrix A0 = new SRMatrix(fsT3.GetNdof());
		double[][] RHS0 = new double[2][fsT3.GetNdof()];
		fsT3.assembleStiff(kxxFunc, param, BVPType.CONSTITUTE_MATRIX_COMMON, A0);
		SRMatrix A = A0.copy();
		fsT3.assembleGlobalVector(kxxFunc, param, SOTS2DT361::elementkxxFirstOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS0);
		double[][] RHS = new double[2][];
		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_Periodic(A, RHS[0], Direct.Y, 1, 3);
		fsT3.applyBC_Periodic(A, RHS[0], Direct.X, 2, 4);
		fsT3.applyBC_Periodic(RHS[1], Direct.Y, 1, 3);
		fsT3.applyBC_Periodic(RHS[1], Direct.X, 2, 4);
		NewIterSSolver solver = new NewIterSSolver(A);
		solver.PCGSSOR(RHS[0], FOCS[0], 1.5, 1);
		solver.PCGSSOR(RHS[1], FOCS[1], 1.5, 1);
		A.clear();

		A = A0.copy();
		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(A, Direct.All, 1, 2, 3, 4);
		solver.setMatrix(A);

		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 2, 3, 4);
		solver.PCGSSOR(RHS[0], FOCS[2], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 2, 3, 4);
		solver.PCGSSOR(RHS[1], FOCS[3], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 2);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 1, 3, 4);
		solver.PCGSSOR(RHS[0], FOCS[4], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 2);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 1, 3, 4);
		solver.PCGSSOR(RHS[1], FOCS[5], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 3);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 1, 2, 4);
		solver.PCGSSOR(RHS[0], FOCS[6], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 3);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 1, 2, 4);
		solver.PCGSSOR(RHS[1], FOCS[7], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 4);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 1, 2, 3);
		solver.PCGSSOR(RHS[0], FOCS[8], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 4);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 1, 2, 3);
		solver.PCGSSOR(RHS[1], FOCS[9], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 3, 4);
		solver.PCGSSOR(RHS[0], FOCS[10], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 3, 4);
		solver.PCGSSOR(RHS[1], FOCS[11], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 3);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 2, 4);
		solver.PCGSSOR(RHS[0], FOCS[12], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 3);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 2, 4);
		solver.PCGSSOR(RHS[1], FOCS[13], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 4);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 2, 3);
		solver.PCGSSOR(RHS[0], FOCS[14], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 4);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 2, 3);
		solver.PCGSSOR(RHS[1], FOCS[15], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 2, 3);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 1, 4);
		solver.PCGSSOR(RHS[0], FOCS[16], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 2, 3);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 1, 4);
		solver.PCGSSOR(RHS[1], FOCS[17], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 2, 4);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 1, 3);
		solver.PCGSSOR(RHS[0], FOCS[18], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 2, 4);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 1, 3);
		solver.PCGSSOR(RHS[1], FOCS[19], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 3, 4);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 1, 2);
		solver.PCGSSOR(RHS[0], FOCS[20], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 3, 4);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 1, 2);
		solver.PCGSSOR(RHS[1], FOCS[21], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 4);
		solver.PCGSSOR(RHS[0], FOCS[22], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 4);
		solver.PCGSSOR(RHS[1], FOCS[23], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 4);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 3);
		solver.PCGSSOR(RHS[0], FOCS[24], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 4);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 3);
		solver.PCGSSOR(RHS[1], FOCS[25], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 3, 4);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 2);
		solver.PCGSSOR(RHS[0], FOCS[26], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 3, 4);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 2);
		solver.PCGSSOR(RHS[1], FOCS[27], 1.5, 1);

		RHS[0] = RHS0[0].clone();
		RHS[1] = RHS0[1].clone();
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 2, 3, 4);
		fsT3.applyBC_MBN_Vector(RHS[0], Direct.All, FOCS[0], 1);
		solver.PCGSSOR(RHS[0], FOCS[28], 1.5, 1);

		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 2, 3, 4);
		fsT3.applyBC_MBN_Vector(RHS[1], Direct.All, FOCS[1], 1);
		solver.PCGSSOR(RHS[1], FOCS[29], 1.5, 1);
		return FOCS;
	}

	public double[][] HomogenizedKxx(MatrixFunc kxxFunc, double[][][] param, double[][] FOCS) {
		double[][] D = new double[2][2];
		double[][] grad = new double[3][2];
		double[][] xy = new double[3][2];
		Mesh2T3 mesh = (Mesh2T3) fsT3.getMesh();
		double[] NN1 = new double[3];
		double[] NN2 = new double[3];
		double N1Dx, N1Dy, N2Dx, N2Dy;
		double wholeArea = 0.0;
		MVO.fill(D, 0.0);
		for (int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, xy);
			MVO.getSlice(FOCS[0], ele, NN1);
			MVO.getSlice(FOCS[1], ele, NN2);
			double area = CG.gradient2DT3(xy, grad);
			double[] dx = MVO.col(grad, 0);
			double[] dy = MVO.col(grad, 1);
			N1Dx = MVO.dot_product(NN1, dx);
			N1Dy = MVO.dot_product(NN1, dy);
			N2Dx = MVO.dot_product(NN2, dx);
			N2Dy = MVO.dot_product(NN2, dy);
			double[] dN1 = {N1Dx, N1Dy};
			double[] dN2 = {N2Dx, N2Dy};
			double[] ip = mesh.getinteriorPointInElement(xy);
			double[][] kxx = kxxFunc.action(ip, mesh.getElementLabel(i), param);
			D[0][0] += (kxx[0][0] + MVO.dot_product(kxx[0], dN1)) * area;
			D[0][1] += (kxx[0][1] + MVO.dot_product(kxx[0], dN2)) * area;
			D[1][0] += (kxx[1][0] + MVO.dot_product(kxx[1], dN1)) * area;
			D[1][1] += (kxx[1][1] + MVO.dot_product(kxx[1], dN2)) * area;
			wholeArea += area;
		}
		MVO.mul(D, 1.0 / wholeArea);
		return D;
	}

	public void KxxSecondOrderCellSolution(SRMatrix A, MatrixFunc kxxFunc, double[][][] param, double[][] FOCS, double[][] D) {
		double[][] SOCS = new double[4][fsT3.GetNdof()];
		double[][] RHS = new double[4][fsT3.GetNdof()];
		DirectSSolver solver = new DirectSSolver(A);
		solver.Factorize_Cholesky(1);
		fsT3.assembleGlobalVector(kxxFunc, param, FOCS, D, SOTS2DT361::elementKxxSecondOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
		MVO.fill(SOCS, 0.0);
		for(int i = 0; i < 4; i++) {
			fsT3.applyBC_MBN(RHS[i], Direct.All, 0.0, 1, 2, 3, 4);
			SOCS[i] = solver.Solve_Cholesky(RHS[i], 1);
		}
	}
	/**
	 * ��������ȴ��������ͷ����ȴ���ϵ����һ�ס����׵������Լ����Ȼ��ȴ���ϵ����
	 * 
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 * @param SOCS
	 */
	public void KxxCellSolutionPeriodic(MatrixFunc kxxFunc, double[][][] param, double[][] FOCS, double[][] D,
			double[][] SOCS) {
		SRMatrix A = new SRMatrix(fsT3.GetNdof());
		fsT3.assembleStiff(kxxFunc, param, BVPType.CONSTITUTE_MATRIX_COMMON, A);
		double[][] RHS = new double[2][fsT3.GetNdof()];
		fsT3.assembleGlobalVector(kxxFunc, param, SOTS2DT361::elementkxxFirstOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
//		System.out.println("max: " + MVO.max(RHS[0]) + "\t " + "min: " + MVO.min(RHS[0]));
//		System.out.println("max: " + MVO.max(RHS[1]) + "\t " + "min: " + MVO.min(RHS[1]));
		//fsT3.applyBC_MBN(A, Direct.All, 1, 2, 3, 4);
		//fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3, 4);
		//fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3, 4);
		//SRMatrix A0 = A.copy();
		//System.out.println(A.display());
		fsT3.applyBC_Periodic(A, RHS[0], Direct.Y, 1, 3);
		fsT3.applyBC_Periodic(A, RHS[0], Direct.X, 2, 4);
		fsT3.applyBC_Periodic(RHS[1], Direct.Y, 1, 3);
		fsT3.applyBC_Periodic(RHS[1], Direct.X, 2, 4);
		MVO.fill(FOCS, 0.0);
//		DirectSSolver solver = new DirectSSolver(A);
		NewIterSSolver solver = new NewIterSSolver(A);
		System.out.println(A.getElement(0, 0));
		solver.PCGSSOR(RHS[0], FOCS[0], 1.5, 1);
		solver.PCGSSOR(RHS[1], FOCS[1], 1.5, 1);
		//solver.CG(RHS[0], FOCS[0], 1);
		//solver.CG(RHS[1], FOCS[1], 1);
//		FOCS[0] = solver.LUSolve(RHS[0], 1, 1e-10);
//		FOCS[1] = solver.LUSolve(RHS[1], 1, 1e-10);
		//solver.PLUSolve(RHS[1], FOCS[1], 1);

		computeHomogenizedKxx(kxxFunc, param, FOCS, D);
		
//		RHS = new double[4][fsT3.GetNdof()];
//		fsT3.assembleGlobalVector(kxxFunc, param, FOCS, D, SOTS2DT361::elementKxxSecondOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
//		MVO.fill(SOCS, 0.0);
//		for(int i = 0; i < 4; i++) {
//			fsT3.applyBC_Periodic(RHS[i], Direct.Y, 1, 3);
//			fsT3.applyBC_Periodic(RHS[i], Direct.X, 2, 4);
//			solver.PCGSSOR(RHS[i], SOCS[i], 1.5, 1);
//		}

		/*
		RHS = new double[4][fsT3.GetNdof()];
		fsT3.assembleGlobalVector(kxxFunc, param, FOCS, D, SOTS2DT361::elementKxxSecondOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
		MVO.fill(SOCS, 0.0);
		fsT3.applyBC_MBN(A0, Direct.All, 1, 2, 3, 4);
		solver.setMatrix(A0);
		for(int i = 0; i < 4; i++) {
			fsT3.applyBC_MBN(RHS[i], Direct.All, 0.0, 1, 2, 3, 4);
			solver.PCGSSSOR(RHS[i], SOCS[i], 1.5, 1);
		} 
		 */
	}
	
	public void EigenValueCellSolution(MatrixFunc kxxFunc, ScalarFunc[] rhoFunc, double[][][] param, double[][] FOCS, double[][] D,
			double[] homoRho, double[][] SOCS) {
		SRMatrix A = new SRMatrix(fsT3.GetNdof());
		fsT3.assembleStiff(kxxFunc, param, BVPType.CONSTITUTE_MATRIX_COMMON, A);
		double[][] RHS = new double[2][fsT3.GetNdof()];
		fsT3.assembleGlobalVector(kxxFunc, param, SOTS2DT361::elementkxxFirstOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
		fsT3.applyBC_MBN(A, Direct.All, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3, 4);
		fsT3.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3, 4);
		MVO.fill(FOCS, 0.0);
		DirectSSolver solver = new DirectSSolver(A);
//		NewIterSSolver solver = new NewIterSSolver(A);
		FOCS[0] = solver.CholeskySolve(RHS[0], 1);
		FOCS[1] = solver.CholeskySolve(RHS[1], 1);

		computeHomogenizedKxx(kxxFunc, param, FOCS, D);
	
		double[] iaValue = fsT3.integrate(rhoFunc, null);
		for(int i = 0, n = iaValue.length; i < n; i++) {
			homoRho[i] = iaValue[i];
		}
		
		RHS = new double[4][fsT3.GetNdof()];
		double[][] DD = new double[3][];
		DD[0] = D[0];
		DD[1] = D[1];
		DD[2] = homoRho;
		fsT3.assembleGlobalVector(kxxFunc, rhoFunc, param, FOCS, DD, SOTS2DT361::elementEigenValueSecondOrder, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
		MVO.fill(SOCS, 0.0);
		for(int i = 0; i < 4; i++) {
			fsT3.applyBC_MBN(RHS[i], Direct.All, 0.0, 1, 2, 3, 4);
			SOCS[i] = solver.CholeskySolve(RHS[i], 1);
		}
	}

	/**
	 * �����ȴ���ϵ��һ�׵��������Ҷ���������
	 * 
	 * @param xy
	 * @param coef
	 * @param tp
	 * @param eleXY
	 */
	public static void elementkxxFirstOrder(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] eleXY) {
		double[][] grad = new double[3][2];
		double area = CG.gradient2DT3(xy, grad);
		//double temp = MVO.sum(coef[0]) * area / 3.0;

		eleXY[0][0] = -MVO.dot_product(coef[0], grad[0]) * area;
		eleXY[0][1] = -MVO.dot_product(coef[0], grad[1]) * area;
		eleXY[0][2] = -MVO.dot_product(coef[0], grad[2]) * area;
		eleXY[1][0] = -MVO.dot_product(coef[1], grad[0]) * area;
		eleXY[1][1] = -MVO.dot_product(coef[1], grad[1]) * area;
		eleXY[1][2] = -MVO.dot_product(coef[1], grad[2]) * area;
		//System.out.println(MVO.toString(coef));
		//System.exit(0);
	}

	/**
	 * �����ȴ���ϵ�����׵����������Ҷ���������
	 * 
	 * @param xy
	 * @param kxx
	 * @param FOCS
	 * @param homo
	 * @param tp
	 * @param eleXY
	 */
	public static void elementKxxSecondOrder(int[] ele, double[][] xy, double[][] kxx, double[][] FOCS, double[][] homo, BVPType tp,
			double[][] eleXY) {
		double[][] grad = new double[3][2];
		double area = CG.gradient2DT3(xy, grad);
		double[] dx = MVO.col(grad, 0);
		double[] dy = MVO.col(grad, 1);
		double N1Dx = MVO.dot_product(dx, FOCS[0]);
		double N1Dy = MVO.dot_product(dy, FOCS[0]);
		double N2Dx = MVO.dot_product(dx, FOCS[1]);
		double N2Dy = MVO.dot_product(dy, FOCS[1]);
		double[] dN1 = {N1Dx, N1Dy};
		double[] dN2 = {N2Dx, N2Dy};
		
		for (int i = 0; i < 3; i++) {
			eleXY[0][i] = (kxx[0][0] - homo[0][0] + MVO.dot_product(kxx[0], dN1) - MVO.dot_product(kxx[0], grad[i]) * MVO.sum(FOCS[0])) * area / 3.0; // 11
			eleXY[1][i] = (kxx[0][1] - homo[0][1] + MVO.dot_product(kxx[0], dN2) - MVO.dot_product(kxx[0], grad[i]) * MVO.sum(FOCS[1])) * area / 3.0; // 12
			eleXY[2][i] = (kxx[1][0] - homo[1][0] + MVO.dot_product(kxx[1], dN1) - MVO.dot_product(kxx[1], grad[i]) * MVO.sum(FOCS[0])) * area / 3.0; // 21
			eleXY[3][i] = (kxx[1][1] - homo[1][1] + MVO.dot_product(kxx[1], dN2) - MVO.dot_product(kxx[1], grad[i]) * MVO.sum(FOCS[1])) * area / 3.0; // 22
		}
//		System.out.println(MVO.toString(xy));
//		System.out.println(Arrays.toString(dx));
//		System.out.println(Arrays.toString(dy));
//		System.out.println(MVO.toString(eleXY));
//		System.exit(0);
	}
	
	public static void elementEigenValueSecondOrder(int[] ele, double[][] xy, double[][] kxx, double[][] rho, double[][] FOCS, double[][] homo, BVPType tp,
			double[][] eleXY) {
		double[][] grad = new double[3][2];
		double area = CG.gradient2DT3(xy, grad);
		double[] dx = MVO.col(grad, 0);
		double[] dy = MVO.col(grad, 1);
		double N1Dx = MVO.dot_product(dx, FOCS[0]);
		double N1Dy = MVO.dot_product(dy, FOCS[0]);
		double N2Dx = MVO.dot_product(dx, FOCS[1]);
		double N2Dy = MVO.dot_product(dy, FOCS[1]);
		double[] dN1 = {N1Dx, N1Dy};
		double[] dN2 = {N2Dx, N2Dy};
		double sumRho = MVO.sum(rho[0]);
		double[] vRho = {0.25 * (sumRho + rho[0][0]), 0.25 * (sumRho + rho[0][1]), 0.25 * (sumRho + rho[0][2])};
		double hRho = homo[2][0];
		for (int i = 0; i < 3; i++) {
			eleXY[0][i] = (kxx[0][0] - vRho[i] * homo[0][0] / hRho + MVO.dot_product(kxx[0], dN1) - MVO.dot_product(kxx[0], grad[i]) * MVO.sum(FOCS[0])) * area / 3.0; // 11
			eleXY[1][i] = (kxx[0][1] - vRho[i] * homo[0][1] / hRho + MVO.dot_product(kxx[0], dN2) - MVO.dot_product(kxx[0], grad[i]) * MVO.sum(FOCS[1])) * area / 3.0; // 12
			eleXY[2][i] = (kxx[1][0] - vRho[i] * homo[1][0] / hRho + MVO.dot_product(kxx[1], dN1) - MVO.dot_product(kxx[1], grad[i]) * MVO.sum(FOCS[0])) * area / 3.0; // 21
			eleXY[3][i] = (kxx[1][1] - vRho[i] * homo[1][1] / hRho + MVO.dot_product(kxx[1], dN2) - MVO.dot_product(kxx[1], grad[i]) * MVO.sum(FOCS[1])) * area / 3.0; // 22
		}
//		System.out.println(MVO.toString(xy));
//		System.out.println(Arrays.toString(dx));
//		System.out.println(Arrays.toString(dy));
//		System.out.println(MVO.toString(eleXY));
//		System.exit(0);
	}

	/**
	 * ��һ�׵������������ȴ������Ȼ�ϵ������
	 * 
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 */
	public void computeHomogenizedKxx(MatrixFunc kxxFunc, double[][][] param, double[][] FOCS, double[][] D) {
		double[][] grad = new double[3][2];
		double[][] xy = new double[3][2];
		Mesh2T3 mesh = (Mesh2T3) fsT3.getMesh();
		double[] NN1 = new double[3];
		double[] NN2 = new double[3];
		double N1Dx, N1Dy, N2Dx, N2Dy;
		double wholeArea = 0.0;
		MVO.fill(D, 0.0);
		for (int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, xy);
			MVO.getSlice(FOCS[0], ele, NN1);
			MVO.getSlice(FOCS[1], ele, NN2);
			double area = CG.gradient2DT3(xy, grad);
			double[] dx = MVO.col(grad, 0);
			double[] dy = MVO.col(grad, 1);
			N1Dx = MVO.dot_product(NN1, dx);
			N1Dy = MVO.dot_product(NN1, dy);
			N2Dx = MVO.dot_product(NN2, dx);
			N2Dy = MVO.dot_product(NN2, dy);
			double[] dN1 = {N1Dx, N1Dy};
			double[] dN2 = {N2Dx, N2Dy};
			double[] ip = mesh.getinteriorPointInElement(xy);
			double[][] kxx = kxxFunc.action(ip, mesh.getElementLabel(i), param);
			D[0][0] += (kxx[0][0] + MVO.dot_product(kxx[0], dN1)) * area;
			D[0][1] += (kxx[0][1] + MVO.dot_product(kxx[0], dN2)) * area;
			D[1][0] += (kxx[1][0] + MVO.dot_product(kxx[1], dN1)) * area;
			D[1][1] += (kxx[1][1] + MVO.dot_product(kxx[1], dN2)) * area;
			wholeArea += area;
		}
		MVO.mul(D, 1.0 / wholeArea);
	}

	public void assembleHeatSOTSSolution(FES2T31 fsFine, double[] T_0, double[] T_1, double[] T_2, FES2T61 fsHomo,
                                         double[] T_Homo, FES2T31 fsCell, double[][] FOCS, double[][] SOCS, int N, double ep) {
		double[] xy = new double[2];
		double[] cellXY = new double[2];
		
		double[] eleVal6 = new double[6];
		double[] eleVal3 = new double[3];

		double[][] grad = new double[2][fsHomo.GetNdof()];
		double[][] Se = new double[3][fsHomo.GetNdof()];
		//fsHomo.computeGradient(T_Homo, grad);
		
		fsHomo.computeFirstAndSecondDerivative(T_Homo, grad, Se);

		double corrector1, corrector2;
		
		DoubleIntegerPair homoData = new DoubleIntegerPair(0.0, 0);
		int[] eleIdx6 = fsHomo.mesh.getElement(0).clone();
		
		DoubleIntegerPair cellData = new DoubleIntegerPair(0.0, 0);
		int[] eleIdx3 = fsCell.mesh.getElement(0).clone();
		double val = 0.0;
		for (int i = 0, n = fsFine.GetNdof(); i < n; i++) {
			fsFine.mesh.getCoordPerNode(i, xy);
			cellXY[0] = xy[0] * N - (int) (xy[0] * N);
			cellXY[1] = xy[1] * N - (int) (xy[1] * N);
			fsHomo.getValue(T_Homo, xy, homoData, eleIdx6, eleVal6);
			T_0[i] = homoData.first;
			fsCell.getValue(FOCS[0], cellXY, cellData, eleIdx3, eleVal3);
			val = cellData.first;
			
//			T_0[i] = fsHomo.getValue(T_Homo, xy, eleIdx6, eleVal6);
//			val = fsCell.getValue(FOCS[0], cellXY, eleIdx3, eleVal3);
			// corrector1 = 0.0;
			// corrector2 = 0.0;
			// corrector1 += val * MVO.dot_product(MVO.getSlice(grad[0], eleIdx6), eleVal6);
			// corrector1 += MVO.dot_product(MVO.getSlice(FOCS[1], eleIdx3), eleVal3)
			// * MVO.dot_product(MVO.getSlice(grad[1], eleIdx6), eleVal6);
			// T_1[i] = T_0[i] + ep * corrector1;
			// corrector2 += MVO.dot_product(MVO.getSlice(SOCS[0], eleIdx3), eleVal3) *
			// MVO.dot_product(MVO.getSlice(Se[0], eleIdx6), eleVal6)
			// + MVO.dot_product(MVO.getSlice(SOCS[1], eleIdx3), eleVal3) *
			// MVO.dot_product(MVO.getSlice(Se[1], eleIdx6), eleVal6)
			// + MVO.dot_product(MVO.getSlice(SOCS[2], eleIdx3), eleVal3) *
			// MVO.dot_product(MVO.getSlice(Se[2], eleIdx6), eleVal6)
			// + MVO.dot_product(MVO.getSlice(SOCS[3], eleIdx3), eleVal3) *
			// MVO.dot_product(MVO.getSlice(Se[3], eleIdx6), eleVal6);

			corrector1 = 0.0;
			corrector1 += val * MVO.dot_product(grad[0], eleIdx6, eleVal6);
			corrector1 += MVO.dot_product(FOCS[1], eleIdx3, eleVal3) * MVO.dot_product(grad[1], eleIdx6, eleVal6);
			T_1[i] = T_0[i] + ep * corrector1;
			
			corrector2 = 0.0;
			corrector2 += MVO.dot_product(SOCS[0], eleIdx3, eleVal3) * MVO.dot_product(Se[0], eleIdx6, eleVal6);
			corrector2 += MVO.dot_product(SOCS[1], eleIdx3, eleVal3) * MVO.dot_product(Se[2], eleIdx6, eleVal6);
			corrector2 += MVO.dot_product(SOCS[2], eleIdx3, eleVal3) * MVO.dot_product(Se[2], eleIdx6, eleVal6);
			corrector2 += MVO.dot_product(SOCS[3], eleIdx3, eleVal3) * MVO.dot_product(Se[1], eleIdx6, eleVal6);

			// Arrays.fill(cellValue, 0.0);
			// Arrays.fill(homoDxxValue, 0.0);
			// for(int j = 0; j < 3; j++) {
			// cellValue[0] += SOCS[0][eleIdx3[j]] * eleVal3[j];
			// cellValue[1] += SOCS[1][eleIdx3[j]] * eleVal3[j];
			// cellValue[2] += SOCS[2][eleIdx3[j]] * eleVal3[j];
			// cellValue[3] += SOCS[3][eleIdx3[j]] * eleVal3[j];
			// }
			// for(int j = 0; j < 6; j++) {
			// homoDxxValue[0] += Se[0][eleIdx6[j]] * eleVal6[j];
			// homoDxxValue[1] += Se[1][eleIdx6[j]] * eleVal6[j];
			// homoDxxValue[2] += Se[2][eleIdx6[j]] * eleVal6[j];
			// homoDxxValue[3] += Se[3][eleIdx6[j]] * eleVal6[j];
			// }
			T_2[i] = T_1[i] + ep * ep * corrector2;
		}
		System.gc();
	}

	public void assembleEigenValueSOTSSolution(FES2T31 fsFine, double[][] u0, double[][] u1, double[][] u2,
                                               ScalarFunc rhoFunc, double[] homoD, double[] d1, double[] d2, double ep) {
		int numEigen = u0.length;
		int nv = fsFine.mesh.getNv();
		double[] rhoVal = fsFine.valueFromFunc(rhoFunc, null);
		double[] u0u0 = new double[nv];
		double[] u0u1 = new double[nv];
		double[] u0u2_u1u1 = new double[nv];
		for(int i = 0; i < numEigen; i++) {
			for(int j = 0; j < nv; j++) {
				u0u0[j] = rhoVal[j] * u0[i][j] * u0[i][j];
				u0u1[j] = rhoVal[j] * (u1[i][j] - u0[i][j]);
				u0u2_u1u1[j] = rhoVal[j] * ((u2[i][j] - u1[i][j]) * u0[i][j] 
	                      + (u1[i][j] - u0[i][j]) * (u1[i][j] - u0[i][j]));
			}
			double lambda1 = - homoD[i] * fsFine.integrate(u0u1) / fsFine.integrate(u0u0) / ep;
			d1[i] = homoD[i] + ep * lambda1;
			d2[i] = d1[i] -homoD[i] * fsFine.integrate(u0u2_u1u1) / fsFine.integrate(u0u0) 
					+ ep * ep * 2 * lambda1 * lambda1 / homoD[i];
		}
	}
	FES2T31 fsT3;
	FES2T61 fsT6;
}
